home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / LinAlg 3.1 / LinAlg / hjmin.cc < prev    next >
Encoding:
Text File  |  1995-12-13  |  6.4 KB  |  215 lines  |  [TEXT/ALFA]
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *              Numerical Math Package
  6.  *        Hook-Jeevse multidimensional minimization
  7.  *
  8.  * Synopsis
  9.  *    double hjmin(b,h,funct)
  10.  *    Vector& b             should be specified to contain the
  11.  *                    initial guess to the point of minimum.
  12.  *                    On return, it contains the location
  13.  *                    of the minimum the function has found
  14.  *    Vector& h            Initial values for steps along
  15.  *                    each direction
  16.  *                    On exit contains the final steps
  17.  *                    before the termination
  18.  *    double f(const Vector& x)    Procedure to compute a function
  19.  *                    value at the specified point
  20.  *
  21.  * The function returns the value of the minimizing function f() at the
  22.  * point of the found minimum.
  23.  *
  24.  * An alternative double hjmin(b,h0,funct) where h0 is a double const
  25.  * uses the same initial steps along each direction. No final steps are
  26.  * reported, though.
  27.  *
  28.  * Algorithm
  29.  *    Hook-Jeevse method of direct search for a function minimum
  30.  *    The method is of the 0. order (i.e. requiring no gradient computation)
  31.  *    See
  32.  *    B.Bondi. Methods of optimization. An Introduction - M.,
  33.  *    "Radio i sviaz", 1988 - 127 p. (in Russian)
  34.  *
  35.  * $Id$
  36.  *
  37.  ************************************************************************
  38.  */
  39.  
  40.  
  41. #include "LinAlg.h"
  42. #include "math_num.h"
  43. #include "std.h"
  44.  
  45. /*
  46.  *------------------------------------------------------------------------
  47.  *        Class to operate on the points in the space (x,f(x))
  48.  */
  49.  
  50. class FPoint
  51. {
  52.   Vector& x;                // Point in the function domain
  53.   double fval;                // Function value at the point
  54.   double (*fproc)(const Vector& x);    // Procedure to compute the function
  55.                     // value
  56.   const bool free_x_on_destructing;    // The flag is needed to free the 
  57.                     // dynamic memory allocated when
  58.                     // x rather than reference to x has
  59.                     // been given to construct FPoint
  60.  
  61. public:
  62.   FPoint(Vector& b, double (*f)(const Vector& x));
  63.   FPoint(const FPoint& fp);
  64.   ~FPoint();
  65.  
  66.   FPoint& operator = (const FPoint& fp);
  67.  
  68.   double f() const            { return fval; }
  69.  
  70.   double fiddle_around(const Vector& h);// Examine the function in the
  71.                     // neighborhood of the current point.
  72.                     // h defines the radius of the region
  73.  
  74.                     // Proceed in direction the function
  75.                     // seems decline
  76.   friend void update_in_direction(FPoint& from, FPoint& to);
  77.  
  78.                     // Decide whether the region embracing
  79.                     // the local min is small enough
  80.   bool is_step_relatively_small(const Vector& h, const double tau);
  81. };
  82.  
  83.                 // Constructor FPoint from array b
  84. inline FPoint::FPoint(Vector& b, double (*f)(const Vector& x))
  85.     : x(b), fproc(f), free_x_on_destructing(false)
  86. {
  87.   fval = (*fproc)(x);
  88. }
  89.  
  90.                 // Constructor by example
  91. FPoint::FPoint(const FPoint& fp)
  92.     : x(*(new Vector(fp.x))), free_x_on_destructing(true)
  93. {
  94.   x = fp.x;
  95.   fval = fp.fval;
  96.   fproc = fp.fproc;
  97. }
  98.                 // Destructor
  99. FPoint::~FPoint()
  100. {
  101.   if( free_x_on_destructing )
  102.     delete &x;
  103. }
  104.  
  105.                 // Assignment; fproc is assumed the same and
  106.                 // is not copied
  107. inline FPoint& FPoint::operator = (const FPoint& fp)
  108. {
  109.   x = fp.x;
  110.   fval = fp.fval;
  111.   return *this;
  112. }
  113.  
  114. /*
  115.  * Examine the function f in the vicinity of the current point
  116.  * by making tentative steps fro/back along each coordinate.
  117.  * Should function decrease, the point is updated to locate the
  118.  * new local min.
  119.  * Examination() returns the minimal function value found in
  120.  * the region.
  121.  *
  122.  */
  123. double FPoint::fiddle_around(const Vector& h)
  124. {
  125.             // Perform a step along a coordinate
  126.   for(register int i=x.q_lwb(); i<=x.q_upb(); i++)
  127.   {
  128.     const double hi = h(i);
  129.     register double xi_old = x(i);      // Old value of x[i]
  130.     register double fnew;
  131.  
  132.     if( x(i) = xi_old + hi, fnew = (*fproc)(x), fnew < fval )
  133.       fval = fnew;            // Step caused f to decrease, OK
  134.     else if( x(i) = xi_old - hi, fnew = (*fproc)(x), fnew < fval )
  135.       fval = fnew;
  136.     else                // No function decline has been
  137.       x(i) = xi_old;                    // found along this coord, back up
  138.   }
  139.   return fval;
  140. }                                                
  141.  
  142.                 // Proceed in the direction the function
  143.                 // seems to fall
  144.                 // to_new = (to - from) + to
  145.                 // from = to (before modification)
  146. void update_in_direction(FPoint& from, FPoint& to)
  147. {
  148.   register int i;
  149.   for(i=(from.x).q_lwb(); i<=(from.x).q_upb(); i++)
  150.   {
  151.     register double t = (to.x)(i);
  152.     (to.x)(i)  += (t - (from.x)(i));
  153.     (from.x)(i) = t;
  154.   }
  155.   from.fval = to.fval;
  156.   to.fval = (*(to.fproc))(to.x);
  157. }
  158.  
  159.                 // Estimate if the point of minimum has
  160.                 // been located accurately enough
  161. bool FPoint::is_step_relatively_small(const Vector& h, const double tau)
  162. {
  163.   register bool it_is_small = true;
  164.   register int i;
  165.   for(i=h.q_lwb(); it_is_small && i<=h.q_upb(); i++)
  166.     it_is_small &= ( h(i) /(1 + fabs(x(i))) < tau );
  167.   return it_is_small;
  168. }
  169.  
  170. /*
  171.  *------------------------------------------------------------------------
  172.  *                Root module
  173.  */
  174.  
  175. double hjmin(Vector& b, Vector& h, double (*ff)(const Vector& x))
  176. {
  177.                 // Function Parameters
  178.   const double tau = 10*EPSILON;        // Termination criterion
  179.   const double threshold = 1e-8;        // Threshold for the function
  180.                                         // decay to be treated as
  181.                                         // significant
  182.   const double step_reduce_factor = 10;   
  183.  
  184.   are_compatible(b,h);
  185.  
  186.   FPoint pmin(b,ff);            // Point of min
  187.   FPoint pbase(pmin);            // Base point
  188.  
  189.   for(;;)            // Main iteration loop
  190.   {                         // pmin is the approximation to min so far
  191.     if( pbase.fiddle_around(h) < pmin.f() - threshold )
  192.     {                         // Function value dropped significantly
  193.       do                                  // from pmin to the point pbase
  194.     update_in_direction(pmin,pbase);// Keep going in the same direction
  195.       while( pbase.fiddle_around(h) < pmin.f() - threshold ); // while it works
  196.       pbase = pmin;            // Save the best approx found
  197.     }
  198.     else                                   // Function didn't fall significantly
  199.       if(                               // upon wandering around pbas
  200.      h *= 1/step_reduce_factor,    // Try to reduce the step then
  201.      pbase.is_step_relatively_small(h,tau) )
  202.     return pmin.f();
  203.   }
  204. }
  205.  
  206.                 // The same as above with the only difference
  207.                 // initial steps are given to be the same
  208.                 // along every direction. The final steps
  209.                 // aren't reported back though
  210. double hjmin(Vector& b,    const double h0, double (*f)(const Vector& x))
  211. {
  212.   Vector h(b.q_lwb(),b.q_upb()); h = h0;
  213.   return hjmin(b,h,f);
  214. }
  215.